library(tidyverse)
library(ggdist)
library(ggside)
library(easystats)
library(patchwork)
library(brms)
library(lavaan)
library(ggraph)
library(tidySEM)
perceptual <- read.csv("../data/scores_perceptual.csv") |>
select(-contains("Sigma")) |>
select(-contains("Beta")) |>
# mutate(across(ends_with("_Error"), \(x) -1*x)) |> # So that higher score = More errors
mutate(across(ends_with("Intercept_Error"), \(x) -1*x)) |> # Not the intercept (re-reverse)
mutate(across(ends_with("_RT"), \(x) -1*x)) |> # So that higher score = Longer RTs
mutate(across(ends_with("Intercept_RT"), \(x) -1*x)) # Not the intercept (re-reverse)
for(i in names(perceptual)[names(perceptual) != "Participant"]) {
perceptual[check_outliers(perceptual[[i]], method = "zscore_robust", threshold = 5), i] <- NA
}
perceptual <- standardize(perceptual) |>
mutate(across(where(is.numeric), as.numeric))
perceptual |>
estimate_density(method = "kernSmooth") |>
separate(Parameter, into = c("Task", "Type", "Parameter", "Model")) |>
ggplot(aes(x=x, y=y, color=Type, linetype = Parameter)) +
geom_line(linetype = 2) +
facet_grid(Parameter~Model)
cor <- correlation::correlation(perceptual, redundant = TRUE, p_adjust = "none")
p_data <- cor |>
mutate(Parameter1 = str_replace_all(str_remove(Parameter1, "Perception_"), "_", " "),
Parameter2 = str_replace_all(str_remove(Parameter2, "Perception_"), "_", " ")) |>
correlation::cor_sort(hclust_method = "ward.D2") |>
cor_lower() |>
mutate(
Text = insight::format_value(r, zap_small = TRUE, digits = 3),
Text = str_replace(str_remove(Text, "^0+"), "^-0+", "-"),
Text = paste0(Text, insight::format_p(p, stars_only = TRUE)),
Parameter2 = fct_rev(Parameter2)
)
p_data |>
ggplot(aes(x = Parameter2, y = Parameter1)) +
geom_tile(aes(fill = r)) +
geom_text(aes(label = Text), size = rel(2), alpha=2/3) +
scale_fill_gradient2(low = "#2196F3", mid = "white", high = "#F44336", midpoint = 0, limit = c(-1, 1), space = "Lab", name = "Correlation", guide = "legend") +
scale_x_discrete(expand = c(0, 0)) +
scale_y_discrete(expand = c(0, 0)) +
labs(title = "Correlation Matrix of Perceptual Task Scores", x = NULL, y = NULL) +
theme_minimal() +
theme(
axis.text.x = element_text(angle = 90, hjust = 1),
plot.title = element_text(hjust = 0.5, face = "bold"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
)
# heatmap(as.matrix(cor))
rez <- parameters::n_factors(select(perceptual, -Participant), n_max=15)
rez
## # Method Agreement Procedure:
##
## The choice of 3 dimensions is supported by 4 (21.05%) methods out of 19 (CNG, Optimal coordinates, Parallel analysis, Kaiser criterion).
plot(rez)
fa <- parameters::factor_analysis(select(perceptual, -Participant),
# cor = as.matrix(cor),
n = 3,
rotation = "oblimin",
fm = "mle",
sort = TRUE
)
insight::print_md(fa)
| Variable | ML2 | ML1 | ML3 | Complexity | Uniqueness |
|---|---|---|---|---|---|
| Perception_Ebbinghaus_Intercept_RT | 1.01 | 0.01 | 0.04 | 1.00 | 4.44e-03 |
| Perception_Ebbinghaus_Difficulty_RT | 1.01 | 5.21e-03 | 0.03 | 1.00 | 5.88e-03 |
| Perception_VerticalHorizontal_Intercept_RT | 0.90 | -0.03 | -0.03 | 1.00 | 0.15 |
| Perception_VerticalHorizontal_Difficulty_RT | 0.89 | -0.05 | -0.03 | 1.01 | 0.16 |
| Perception_MullerLyer_Intercept_RT | 0.89 | 0.03 | -0.01 | 1.00 | 0.22 |
| Perception_MullerLyer_Difficulty_RT | 0.88 | 0.02 | -0.02 | 1.00 | 0.23 |
| Perception_VerticalHorizontal_Difficulty_Error | 2.54e-03 | 1.00 | -9.37e-03 | 1.00 | 2.71e-03 |
| Perception_VerticalHorizontal_Intercept_Error | 1.84e-03 | 1.00 | -8.15e-03 | 1.00 | 2.70e-03 |
| Perception_Ebbinghaus_Difficulty_Error | -0.08 | 0.41 | 0.24 | 1.70 | 0.64 |
| Perception_Ebbinghaus_Intercept_Error | -0.05 | 0.36 | 0.24 | 1.78 | 0.71 |
| Perception_MullerLyer_Intercept_Error | 8.08e-03 | -1.60e-03 | 1.00 | 1.00 | 0.02 |
| Perception_MullerLyer_Difficulty_Error | -8.86e-03 | 1.47e-04 | 0.99 | 1.00 | 5.00e-03 |
The 3 latent factors (oblimin rotation) accounted for 82.17% of the total variance of the original data (ML2 = 43.72%, ML1 = 20.12%, ML3 = 18.33%).
plot(fa) +
theme_minimal()
# psych::omega(data, fm = "mle", nfactors=7)
# fa <- psych::fa.multi(as.matrix(cor), nfactors = 3, nfact2 = 1, n.obs = nrow(random))
# psych::fa.multi.diagram(fa)
data <- predict(fa)
# rez <- parameters::n_factors(data, rotation = "varimax")
# plot(rez)
fa2 <- parameters::factor_analysis(data,
n = 1,
rotation = "varimax",
fm = "mle",
sort = TRUE
)
insight::print_md(fa2)
| Variable | ML1 | Complexity | Uniqueness |
|---|---|---|---|
| ML3 | 0.77 | 1.00 | 0.41 |
| ML1 | 0.66 | 1.00 | 0.57 |
| ML2 | -0.42 | 1.00 | 0.82 |
The unique latent factor (varimax rotation) accounted for 40.10% of the total variance of the original data.
fit_sem <- function(model, data, ...) {
sem(model = model,
data = data[complete.cases(data), ],
bounds = "standard",
estimator = "ML",
auto.fix.first=TRUE,
control=list(iter.max=200000, eval.max = 20000),
...)
}
names <- names(select(perceptual, -Participant))
data <- setNames(as.data.frame(str_split_fixed(names, "_", 4)), c("Task", "Illusion_Type", "Parameter", "Outcome"))
data$Name <- names
# Model 0 - by EFA
by_efa <- efa_to_cfa(fa, threshold = "max")
# cat(by_model)
m0 <- fit_sem(by_efa, data=perceptual)
# Model 1 - by Model
by_model <- paste(sort(paste0("p_", ifelse(str_detect(data$Outcome, "RT"), "RT", "Error"), " =~ ", data$Name)), collapse = "\n")
# cat(by_model)
m1 <- fit_sem(by_model, data=perceptual)
# Model 2 - by Component
by_component <- paste(sort(paste0("p_", data$Outcome, " =~ ", data$Name)), collapse = "\n")
# cat(by_component)
m2 <- fit_sem(by_component, data=perceptual)
# Model 3 - by Type
by_type <- paste(paste(paste0(data$Illusion_Type), "=~", data$Name), collapse = "\n")
# cat(by_type)
m3 <- fit_sem(by_type, data=perceptual)
# Model 4 - by Parameter
by_param <- paste(paste(paste0(data$Parameter), "=~", data$Name), collapse = "\n")
# cat(by_param)
m4 <- fit_sem(by_param, data=perceptual)
# Model 4 - Unique
m5 <- fit_sem(paste0("p =~ ", paste(names, collapse = " + ")), data=perceptual)
anova(m2, m0, m1, m3, m4, m5)
## Chi-Squared Difference Test
##
## Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
## m0 51 7574 7687 6645
## m2 53 6717 6822 5792 -853 2 1
## m1 53 6717 6822 5792 0 0
## m4 53 10050 10154 9124 3332 0
## m5 54 10273 10373 9349 225 1 <2e-16 ***
## m3 9610 9722
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
test_bf(m2, m0, m1, m3, m4, m5)
## Bayes Factors for Model Comparison
##
## Model BF
## [2] m0 1.28e-188
## [3] m1 1.00
## [4] m3 0.00e+00
## [5] m4 0.00e+00
## [6] m5 0.00e+00
##
## * Against Denominator: [1] m2
## * Bayes Factor Type: BIC approximation
# One factor
m2b <- fit_sem(paste0(by_component, "\np =~ p_RT + p_Error"), data=perceptual)
anova(m2, m2b)
## Chi-Squared Difference Test
##
## Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
## m2b 52 6749 6858 5822
## m2 53 6717 6822 5792 -30.2 1 1
test_bf(m2, m2b)
## Bayes Factors for Model Comparison
##
## Model BF
## [2] m2b 1.28e-08
##
## * Against Denominator: [1] m2
## * Bayes Factor Type: BIC approximation
compare_performance(m2, m2b, metrics = c("RMSEA", "NFI", "CFI"))
## # Comparison of Model Performance Indices
##
## Name | Model | RMSEA | NFI | CFI
## -------------------------------------
## m2 | lavaan | 0.474 | 0.620 | 0.622
## m2b | lavaan | 0.480 | 0.618 | 0.620
summary(m2, standardize = TRUE, fit.measures = TRUE)
## lavaan 0.6-12 ended normally after 887 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 25
##
## Number of observations 481
##
## Model Test User Model:
##
## Test statistic 5791.869
## Degrees of freedom 53
## P-value (Chi-square) 0.000
##
## Model Test Baseline Model:
##
## Test statistic 15236.653
## Degrees of freedom 66
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.622
## Tucker-Lewis Index (TLI) 0.529
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -3333.612
## Loglikelihood unrestricted model (H1) -437.678
##
## Akaike (AIC) 6717.224
## Bayesian (BIC) 6821.621
## Sample-size adjusted Bayesian (BIC) 6742.274
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.474
## 90 Percent confidence interval - lower 0.464
## 90 Percent confidence interval - upper 0.485
## P-value RMSEA <= 0.05 0.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.152
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## p_Error =~
## Prcptn_Ebb_D_E 1.000 0.495 0.520
## Prcptn_Ebb_I_E 1.004 0.110 9.138 0.000 0.497 0.496
## Prcptn_MlL_D_E 0.978 0.107 9.163 0.000 0.484 0.498
## Prcptn_MlL_I_E 0.971 0.107 9.114 0.000 0.481 0.494
## Prcptn_VrH_D_E 1.972 0.148 13.346 0.000 0.976 1.000
## Prcptn_VrH_I_E 1.972 0.148 13.346 0.000 0.976 1.000
## p_RT =~
## Prcptn_Eb_D_RT 1.000 0.983 0.996
## Prcptn_Eb_I_RT 1.003 0.005 215.961 0.000 0.986 0.999
## Prcptn_ML_D_RT 0.858 0.021 40.338 0.000 0.844 0.881
## Prcptn_ML_I_RT 0.866 0.021 41.827 0.000 0.851 0.889
## Prcptn_VH_D_RT 0.925 0.019 48.505 0.000 0.909 0.914
## Prcptn_VH_I_RT 0.929 0.019 49.606 0.000 0.913 0.918
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## p_Error ~~
## p_RT -0.120 0.025 -4.887 0.000 -0.247 -0.247
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .Prc_E_D_E 0.661 0.043 15.508 0.000 0.661 0.730
## .Prc_E_I_E 0.755 0.049 15.508 0.000 0.755 0.754
## .Pr_ML_D_E 0.710 0.046 15.508 0.000 0.710 0.752
## .Pr_ML_I_E 0.714 0.046 15.508 0.000 0.714 0.756
## .Pr_VH_D_E 0.001 0.001 0.756 0.450 0.001 0.001
## .Pr_VH_I_E (lb) 0.000 0.001 0.000 0.000 0.000
## .Pr_E_D_RT 0.008 0.001 6.827 0.000 0.008 0.008
## .Pr_E_I_RT 0.002 0.001 1.531 0.126 0.002 0.002
## .P_ML_D_RT 0.204 0.013 15.425 0.000 0.204 0.223
## .P_ML_I_RT 0.193 0.013 15.418 0.000 0.193 0.210
## .P_VH_D_RT 0.162 0.011 15.380 0.000 0.162 0.164
## .P_VH_I_RT 0.156 0.010 15.373 0.000 0.156 0.158
## p_Error 0.245 0.040 6.130 0.000 1.000 1.000
## p_RT 0.966 0.063 15.381 0.000 1.000 1.000
graph_sem(m2, layout = get_layout(m2, layout_algorithm = "layout_with_lgl"))
scores_p <- mutate(as.data.frame(predict(m2)), Participant = perceptual$Participant[complete.cases(perceptual)])
empirical <- read.csv("../data/preprocessed_perceptual.csv")
empirical <- empirical |>
group_by(Participant, Illusion_Type) |>
summarize(Empirical_Errors = as.numeric(sum(Error, na.rm=TRUE))) |>
full_join(
empirical |>
filter(Error == 0) |>
group_by(Participant, Illusion_Type) |>
summarize(Empirical_MeanRT = mean(RT, na.rm=TRUE)), by = c("Participant", "Illusion_Type")) |>
ungroup() |>
pivot_wider(names_from = "Illusion_Type", values_from = c("Empirical_Errors", "Empirical_MeanRT")) |>
left_join(scores_p, by = "Participant") |>
left_join(perceptual, by = "Participant")
correlation::correlation(select(empirical, starts_with("Empirical")), select(empirical, -Participant, -starts_with("Empirical"))) |>
summary() |>
plot() +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
illusion <- read.csv("../data/scores_illusion1.csv") |>
select(-contains("Diff")) |>
select(-contains("Interaction")) |>
select(-contains("Sigma")) |>
select(-contains("Beta"))
# mutate(across(ends_with("_Error"), \(x) -1*x)) # So that higher score = More errors
for(i in names(illusion)[names(illusion) != "Participant"]) {
illusion[check_outliers(illusion[[i]], method = "zscore_robust", threshold = 5), i] <- NA
}
illusion <- standardize(illusion) |>
mutate(across(where(is.numeric), as.numeric))
illusion |>
estimate_density(method = "kernSmooth") |>
separate(Parameter, into = c("Type", "Parameter", "Model")) |>
ggplot(aes(x=x, y=y, color=Type)) +
geom_line(linewidth=1) +
facet_grid(Model ~ Parameter, scales = "free")
cor <- correlation::correlation(illusion, redundant = TRUE, p_adjust = "none")
p_data <- cor |>
correlation::cor_sort(hclust_method = "ward.D2") |>
cor_lower() |>
mutate(
Text = insight::format_value(r, zap_small = TRUE, digits = 3),
Text = str_replace(str_remove(Text, "^0+"), "^-0+", "-"),
Text = paste0(Text, insight::format_p(p, stars_only = TRUE)),
Parameter2 = fct_rev(Parameter2)
)
p_data |>
ggplot(aes(x = Parameter2, y = Parameter1)) +
geom_tile(aes(fill = r)) +
geom_text(aes(label = Text), size = rel(4), alpha=3/3) +
scale_fill_gradient2(low = "#2196F3", mid = "white", high = "#F44336", midpoint = 0, limit = c(-1, 1), space = "Lab", name = "Correlation", guide = "legend") +
scale_x_discrete(expand = c(0, 0)) +
scale_y_discrete(expand = c(0, 0)) +
labs(title = "Correlation Matrix", x = NULL, y = NULL) +
theme_minimal() +
theme(
axis.text.x = element_text(angle = 90, hjust = 1),
plot.title = element_text(hjust = 0.5, face = "bold"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
)
# heatmap(as.matrix(cor))
rez <- parameters::n_factors(select(illusion, -Participant), n_max=9)
rez
## # Method Agreement Procedure:
##
## The choice of 2 dimensions is supported by 7 (43.75%) methods out of 16 (Optimal coordinates, Parallel analysis, Kaiser criterion, VSS complexity 1, VSS complexity 2, BIC, BIC (adjusted)).
plot(rez)
fa <- parameters::factor_analysis(select(illusion, -Participant),
# cor = as.matrix(cor),
n = 2,
rotation = "oblimin",
fm = "mle",
sort = TRUE
)
insight::print_md(fa)
| Variable | ML2 | ML1 | Complexity | Uniqueness |
|---|---|---|---|---|
| MullerLyer_Strength_RT | 0.88 | -0.04 | 1.00 | 0.25 |
| Ebbinghaus_Strength_RT | 0.77 | -0.12 | 1.04 | 0.43 |
| VerticalHorizontal_Strength_RT | 0.72 | 0.19 | 1.14 | 0.39 |
| Ebbinghaus_Strength_Error | 0.45 | 0.12 | 1.15 | 0.75 |
| MullerLyer_Strength_Error | 0.29 | 0.21 | 1.84 | 0.84 |
| VerticalHorizontal_Strength_Error | 1.95e-04 | 1.00 | 1.00 | 5.00e-03 |
The 2 latent factors (oblimin rotation) accounted for 55.57% of the total variance of the original data (ML2 = 36.66%, ML1 = 18.91%).
plot(fa) +
theme_minimal()
# psych::omega(data, fm = "mle", nfactors=7)
# fa <- psych::fa.multi(as.matrix(cor), nfactors = 3, nfact2 = 1, n.obs = nrow(random))
# psych::fa.multi.diagram(fa)
# data <- predict(fa)
#
# rez <- parameters::n_factors(na.omit(data), rotation = "oblimin")
# plot(rez)
#
# fa2 <- parameters::factor_analysis(data,
# n = 2,
# rotation = "oblimin",
# fm = "mle",
# sort = TRUE
# )
#
# insight::print_md(fa2)
names <- names(select(illusion, -Participant))
data <- setNames(as.data.frame(str_split_fixed(names, "_", 3)), c("Illusion_Type", "Parameter", "Outcome"))
data$Name <- names
# Model 1 - EFA
by_efa <- efa_to_cfa(fa, threshold = "max")
m1 <- fit_sem(by_efa, data=illusion)
# Model 2 - by Model
by_model <- paste(sort(paste(data$Outcome, "=~", data$Name)), collapse = "\n")
# cat(by_model)
m2 <- fit_sem(by_model, data=illusion)
# Model 3 - by Type
by_type <- paste(paste(paste0(data$Illusion_Type), "=~", data$Name), collapse = "\n")
# cat(by_type)
m3 <- fit_sem(by_type, data=illusion)
# Model 4 - Unique
m4 <- fit_sem(paste0("p =~ ", paste(names, collapse = " + ")), data=illusion)
anova(m1, m2, m3, m4)
## Chi-Squared Difference Test
##
## Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
## m3 6 7580 7643 64.3
## m2 8 7598 7653 86.9 22.6 2 1.3e-05 ***
## m1 9 7611 7661 101.4 14.5 1 0.00014 ***
## m4 9 7611 7661 101.4 0.0 0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
test_bf(m1, m2, m3, m4)
## Bayes Factors for Model Comparison
##
## Model BF
## [2] m2 63.13
## [3] m3 1.02e+04
## [4] m4 1.000
##
## * Against Denominator: [1] m1
## * Bayes Factor Type: BIC approximation
m5 <- fit_sem(paste0(by_type, "\ni =~ Ebbinghaus + MullerLyer + VerticalHorizontal"), data=illusion)
anova(m1, m5)
## Chi-Squared Difference Test
##
## Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
## m5 6 7580 7643 64.3
## m1 9 7611 7661 101.4 37.1 3 4.5e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
test_bf(m1, m5)
## Bayes Factors for Model Comparison
##
## Model BF
## [2] m5 1.02e+04
##
## * Against Denominator: [1] m1
## * Bayes Factor Type: BIC approximation
compare_performance(m1, m5, metrics = c("RMSEA", "NFI", "CFI"))
## # Comparison of Model Performance Indices
##
## Name | Model | RMSEA | NFI | CFI
## -------------------------------------
## m1 | lavaan | 0.144 | 0.879 | 0.887
## m5 | lavaan | 0.140 | 0.923 | 0.929
display(performance(m5))
| Chi2(6) | p (Chi2) | Baseline(15) | p (Baseline) | GFI | AGFI | NFI | NNFI | CFI | RMSEA | RMSEA CI | p (RMSEA) | RMR | SRMR | RFI | PNFI | IFI | RNI | Loglikelihood | AIC | BIC | BIC_adjusted |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 64.30 | < .001 | 835.82 | < .001 | 0.96 | 0.85 | 0.92 | 0.82 | 0.93 | 0.14 | [0.11, 0.17] | < .001 | 0.06 | 0.06 | 0.81 | 0.37 | 0.93 | 0.93 | -3774.78 | 7579.57 | 7642.58 | 7594.97 |
summary(m5, standardize = TRUE, fit.measures = TRUE)
## lavaan 0.6-12 ended normally after 78 iterations
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 15
##
## Number of observations 493
##
## Model Test User Model:
##
## Test statistic 64.303
## Degrees of freedom 6
## P-value (Chi-square) 0.000
##
## Model Test Baseline Model:
##
## Test statistic 835.820
## Degrees of freedom 15
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 0.929
## Tucker-Lewis Index (TLI) 0.822
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -3774.784
## Loglikelihood unrestricted model (H1) -3742.632
##
## Akaike (AIC) 7579.568
## Bayesian (BIC) 7642.575
## Sample-size adjusted Bayesian (BIC) 7594.965
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.140
## 90 Percent confidence interval - lower 0.111
## 90 Percent confidence interval - upper 0.172
## P-value RMSEA <= 0.05 0.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.063
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## Ebbinghaus =~
## Ebbnghs_Strn_E 1.000 0.498 0.498
## Ebbnghs_Str_RT 1.493 0.150 9.948 0.000 0.744 0.766
## MullerLyer =~
## MllrLyr_Strn_E 1.000 0.323 0.325
## MllrLyr_Str_RT 2.776 0.457 6.076 0.000 0.898 0.915
## VerticalHorizontal =~
## VrtclHrznt_S_E 1.000 0.340 0.340
## VrtclHrzn_S_RT 2.881 0.530 5.433 0.000 0.980 1.000
## i =~
## Ebbinghaus 1.000 0.963 0.963
## MullerLyer 0.635 0.121 5.250 0.000 0.942 0.942
## VerticalHrzntl 0.532 0.112 4.738 0.000 0.750 0.750
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .Ebbng_S_E 0.752 0.052 14.429 0.000 0.752 0.752
## .Ebbn_S_RT 0.388 0.052 7.466 0.000 0.388 0.413
## .MllrL_S_E 0.885 0.058 15.365 0.000 0.885 0.894
## .MllL_S_RT 0.156 0.092 1.699 0.089 0.156 0.162
## .VrtcH_S_E 0.885 0.058 15.125 0.000 0.885 0.884
## .VrtH_S_RT (lb) 0.000 0.130 0.000 0.000 0.000
## .Ebbinghas 0.018 0.021 0.859 0.390 0.073 0.073
## .MullerLyr 0.012 0.011 1.070 0.285 0.113 0.113
## .VrtclHrzn 0.051 0.014 3.752 0.000 0.438 0.438
## i 0.230 0.045 5.100 0.000 1.000 1.000
# graph_sem(m5, layout = get_layout(m5, layout_algorithm = "layout_with_lgl"))
colors <- c("i" = "#9C27B0",
"Ebbinghaus"="#3F51B5",
"VerticalHorizontal"="#FF5722",
"MullerLyer" = "#4CAF50",
"black"="black")
edges <- tidySEM::get_edges(m5) |>
mutate(color = ifelse(sign(as.numeric(est_std)) >= 0, "1", "-1"),
width = abs(as.numeric(est_std)),
type = ifelse(from == "i", "curved", "straight"),
text = paste(est_std, confint_std)) |>
filter(lhs != rhs)
nodes <- tidySEM::get_nodes(m5) |>
mutate(angle = ifelse(grepl("_", name), 90, 0),
hjust = case_when(
grepl("_", name) ~ 1.5,
name == "i" ~ 0.5,
TRUE ~ 0.5),
size = case_when(
grepl("_", name) ~ 1,
name == "i" ~ 7,
TRUE ~ 6),
textsize = case_when(
grepl("_", name) ~ 1,
name == "i" ~ 3,
TRUE ~ 1.1),
face = case_when(
grepl("_", name) ~ "italic",
name == "i" ~ "bold.italic",
TRUE ~ "plain"),
color = str_remove_all(label, " - .*"),
textcolor = ifelse(grepl("_", name), "black", "white"),
label = case_when(
label == "MullerLyer" ~ "Müller-Lyer",
label == "VerticalHorizontal" ~ "Vertical-\nHorizontal",
str_detect(label, "_Error") ~ "Error",
str_detect(label, "_RT") ~ "RT ",
TRUE ~ label
),
family = ifelse(name == "i", "serif", "sans"))
p_sem <- tidygraph::tbl_graph(nodes = nodes, edges = edges) |>
ggraph(layout = "sugiyama") +
geom_edge_arc(aes(filter=type!="straight",
edge_width = width,
color = color,
label=est_std),
strength = 0.03 * c(-1, 1, 1),
angle_calc="along",
label_dodge=unit(-0.017, "npc"),
label_size = rel(8)) +
geom_edge_link(aes(filter=type=="straight",
edge_width = width,
color = color,
label=est_std),
angle_calc="along",
label_dodge=unit(-0.017, "npc"),
label_size = rel(8)) +
geom_node_point(aes(shape = shape, color=color, size=size)) +
geom_node_text(aes(label = label, angle = angle, hjust=hjust),
color=nodes$textcolor,
fontface=nodes$face,
size = rel(nodes$textsize*10),
family=nodes$family) +
scale_y_continuous(expand = expansion(add=c(0.5, 0.5))) +
scale_edge_width_continuous(range=c(0.3, 6)) +
scale_edge_color_manual(values=c("1"="#263238", "-1"="#B71C1C")) +
scale_color_manual(values=colors, guide="none") +
scale_shape_manual(values=c("oval"="circle", "rect"="square"), guide="none") +
scale_size_continuous(range = c(10, 95), guide="none") +
guides(edge_colour = "none", edge_width = "none") +
labs(title = "Structural Equation Model") +
theme_graph() +
theme(plot.title = element_text(hjust=0.5))
p_sem
ggsave("figures/figure_sem.png", p_sem, width=15, height=15)
# save(p_sem, file="models/p_sem.Rdata")
scores_i <- mutate(as.data.frame(predict(m5)), Participant = illusion$Participant[complete.cases(illusion)])
empirical <- read.csv("../data/preprocessed_illusion1.csv")
empirical <- empirical |>
filter(Illusion_Effect == "Incongruent") |>
group_by(Participant, Illusion_Type) |>
summarize(Empirical_Errors = as.numeric(sum(Error, na.rm=TRUE))) |>
full_join(
empirical |>
filter(Illusion_Effect == "Incongruent", Error == 0) |>
group_by(Participant, Illusion_Type) |>
summarize(Empirical_MeanRT = mean(RT, na.rm=TRUE)), by = c("Participant", "Illusion_Type")) |>
ungroup() |>
pivot_wider(names_from = "Illusion_Type", values_from = c("Empirical_Errors", "Empirical_MeanRT")) |>
left_join(scores_i, by = "Participant") |>
left_join(illusion, by = "Participant")
correlation::correlation(select(empirical, starts_with("Empirical")), select(empirical, -Participant, -starts_with("Empirical"))) |>
summary() |>
plot() +
theme(axis.text.x = element_text(angle = 90, hjust = 1))
write.csv(full_join(left_join(scores_i, read.csv("../data/scores_illusion1.csv"), by = "Participant"), scores_p, by = "Participant"), "../data/scores_sem.csv", row.names = FALSE)